/*
 * MAtoms.cpp
 *
 *  Created on: Aug 7, 2013
 *      Author: marchi
 */

#include "MAtoms.h"



template <typename T, typename U>
class create_map{
private:
    map<T, U> m_map;
public:
    create_map(const T& key, const U& val){m_map[key] = val;}

    create_map<T, U>& operator()(const T& key, const U& val) {
    	m_map[key] = val;
    	return *this;
    }

    operator map<T, U>(){return m_map;}
};

template <typename T>
class create_vector{
private:
    vector<T> m_vec;
public:
    create_vector(const T& data){m_vec.push_back(data);}

    create_vector<T>& operator()(const T& data) {
    	m_vec.push_back(data);
    	return *this;
    }

    operator vector<T>(){return m_vec;}
};

int MAtoms::calls=0;
vector<string> * MAtoms::ResList0=NULL;
vector<int> * MAtoms::ResIndx0=NULL;

map<string,double> MAtoms::MMass=create_map<string,double>
("H",1.008)("C",12.011)("N",14.0067)("O",15.9994)("S",32.06)("P",30.974)
("FE",55.847)("CA",40.080)("ZN",65.37)("NA",22.98977)("MG",24.305);
map<string,double> MAtoms::MMassNCH=create_map<string,double>
("H",0.0)("C",0.0)("N",14.0067)("O",15.9994)("S",32.06)("P",30.974)
("FE",55.847)("CA",40.080)("ZN",65.37)("NA",22.98977)("MG",24.305);
string MAtoms::MIons="CA NA SOD CAL ZN";
string MAtoms::Mdetg="AOT SDS LEA LDA DOP";

class PBCvect{
	static vector<Dvect> * nxyz;
	const int M;
public:
	PBCvect(): M(3){
		nxyz=new vector<Dvect>;
		for(int o=0;o<M;o++){
			double nx=static_cast<double>(o-1);
			for(int p=0;p<M;p++){
				double ny=static_cast<double>(p-1);
				for(int q=0;q<M;q++){
					double nz=static_cast<double>(q-1);
					nxyz->push_back(Dvect(nx,ny,nz));
				}
			}
		}
	}
	static vector<Dvect> & getVec(){return *nxyz;}
};
struct op_iComp {
  bool operator() (double i,double j) { return (i>j);}
};

struct op_sortResPDB{
	bool operator()(string x, string y){
		stringstream ssx(x.substr(21,5));
		stringstream ssy(y.substr(21,5));
		int nx, ny;
		ssx>>nx; ssy>>ny;
		return (nx < ny);};
};
struct op_sortAtnoPDB{
	bool operator()(string x, string y){
		stringstream ss1x(x.substr(6,5));
		stringstream ss1y(y.substr(6,5));
		int nx, ny;
		ss1x>>nx; ss1y>>ny;
		return (nx < ny);
	}
};


class Jacob{
	double ** iner;
	double ** ei;
	double * Im;
	static Dvect Im_s;
	void  eigsort(Dvect  & d, Matrix  * v=NULL){
		int k;
		int n=DIM;
		for (int i=0;i<n-1;i++) {
			double p=d[k=i];
			for (int j=i;j<n;j++)
				if (d[j] >= p) p=d[k=j];
			if (k != i) {
				d[k]=d[i];
				d[i]=p;
				if (v != NULL){
					for (int j=0;j<n;j++) {
						p=(*v)[j][i];
						(*v)[j][i]=(*v)[j][k];
						(*v)[j][k]=p;
					}
				}
			}
		}
	}
public:
	Jacob():  iner(NULL), ei(NULL),Im(NULL){};
	~Jacob(){
		for(int o=0;o<DIM;o++) delete [] iner[o];
		delete [] iner;
		for(int o=0;o<DIM;o++) delete [] ei[o];
		delete [] ei;
		delete [] Im;
	}
	void operator()(Matrix & iner0, Dvect & Im0, Matrix & ei0, int & nrot){
		if(!iner) {
			iner=new double * [DIM];
			for(int o=0;o<DIM;o++) iner[o]=new double [DIM];
		}
		if(!ei) {
			ei=new double * [DIM];
			for(int o=0;o<DIM;o++) ei[o]=new double [DIM];
		}
		if(!Im) Im=new double [DIM];
		for(int o=0;o<DIM;o++){
			for(int p=0;p<DIM;p++){
				iner[o][p]=iner0[o][p];
			}
		}
		jacobi(iner,DIM,Im,ei,&nrot);

		for(int o=0;o<DIM;o++){
			Im0[o]=Im[o];
			for(int p=0;p<DIM;p++){
				ei0[o][p]=ei[o][p];
			}
		}

		eigsort(Im0,&ei0);
	}

};
Dvect Jacob::Im_s=0.0;

vector<Dvect> * PBCvect::nxyz=NULL;

class CellComp{
	Dvect Xd;
	Matrix co;
public:
	CellComp(Dvect& x, Dvect & y, Matrix & c): Xd(x-y), co(c){};
	CellComp(Dvect& x, Matrix & c): Xd(x), co(c){};
	bool operator()(const Dvect & x, const Dvect & y){
		Dvect x1=Xd+x;
		Dvect x2=Xd+y;
		Dvect xc1=co*x1;
		Dvect xc2=co*x2;
		return xc1.Norm() < xc2.Norm();
	}
};
void MAtoms::Percolate(vector<vector<int> > & Mols){
	vector<string> resn=vector<string>(PDB.size());
	for(size_t o=0;o<resn.size();o++){resn[o]=PDB[o].resn;}
	vector<string> atmss=vector<string>(nr);
	for(size_t o=0;o<PDB.size();o++)
		atmss[o]=PDB[o].atn;
	Perco=new Percolation(SelRes,rd,resn,atmss);

}
Percolation & MAtoms::Percolate() {
	try{
		if(!Perco) throw string("Should initialize percolation. Abort.");
	} catch(const string & s){
		cout << s <<endl;
		exit(1);
	}
	vector<Dvect> v(x,x+nr);
	Matrix co(Mt.getCO());
	Matrix oc(Mt.getOC());
	Perco->doContacts(v,co,oc);
	Perco->gCluster();
	Perco->Accumulate();
	return *Perco;
}
MAtoms::MAtoms(): iCountRot{0},Rg_count(0), Perco(NULL), R_cmx{new CenterMass}{};
MAtoms::MAtoms(int nn): iCountRot{0},Rg_count(0), Perco(NULL){
	switch(nn){
	case 1:
		R_cmx=new CenterMassBW;
		break;
	case 2:
		R_cmx=new CenterMassBW2;
		break;
	case 3:
		R_cmx=new CenterMassBW3;
		break;
	default:
		R_cmx=new CenterMass;
	}
}
void MAtoms::initLists(TopolPDB & y,vector<string> & L){
	PDB=vector<PDBdata>(y.Size());

	for(size_t o=0;o<PDB.size();o++) PDB[o]=y[o];
	ResList0=new vector<string>(L);
	ResIndx0=new vector<int>;
	int ia=0;
	for(size_t o=0;o<PDB.size();o++){
		string sub1=PDB[o].atn;
		string sub2=PDB[o].resn;
		atres.push_back(sub2);
		atname.push_back(sub1);
		size_t mst=sub1.find_first_not_of("1234")?1:0;

		if(MIons.find(sub2) == string::npos) {
			mass.push_back(MMass[sub1.substr(mst,1)]);
			massNCH.push_back(MMassNCH[sub1.substr(mst,1)]);
		}
		else {
			mass.push_back(MMass[sub1.substr(mst,2)]);
			massNCH.push_back(MMassNCH[sub1.substr(mst,2)]);
		}
		if(!ResList0 || find(ResList0->begin(),ResList0->end(),sub2) != ResList0->end()) {
			ResIndx0->push_back(ia);
		}
		ia++;
	}
}
void MAtoms::CalcGyro(vector<double> & massa,vector<Gyration> & Rg){
	vector<vector<int> > mCluster=Perco->getCluster();
	vector<vector<int> > mAtoms=Perco->getAtoms();
	for(size_t o=0;o<mCluster.size();o++){
		Rg[o]=0.0;
		Dvect cm=0.0,cmG=0.0;
		double unit_nmm2=1.0/(unit_nm*unit_nm);

		int ntot=0;
		double tmass=0.0;
		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];
			for(size_t q=0;q<mAtoms[n].size();q++){
				int i=mAtoms[n][q];
				for(int o1=0;o1<DIM;o1++) cm[o1]+=x[i][o1]*massa[i];
				for(int o1=0;o1<DIM;o1++) cmG[o1]+=x[i][o1];
				tmass+=massa[i];
				ntot++;
			}
		}
		cm/=tmass;
		cmG/=static_cast<int> (ntot);
		Dvect x0;
		Dvect Gm,Im,axis;
		Matrix Giner=0.0;
		Matrix ei,ei2=0.0;
		double MyRg=0;
		int oa=0;
		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];
			for(size_t q=0;q<mAtoms[n].size();q++){
				int nn=mAtoms[n][q];
				for(int o1=0;o1<DIM;o1++) x0[o1]=cmG[o1]-x[nn][o1];
				Giner+=x0%x0*unit_nmm2;
			}
		}
		Giner/=static_cast<double> (ntot);
		int nrot;
		Jacob()(Giner,Gm,ei,nrot);
		for(int o1=0;o1<DIM;o1++) MyRg+=Gm[o1];
		Matrix Inertia=0.0;
		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];
			for(size_t q=0;q<mAtoms[n].size();q++){
				int nn=mAtoms[n][q];
				for(int o1=0;o1<DIM;o1++) x0[o1]=cm[o1]-x[nn][o1];

				Inertia[XX][XX]+=(sqr(x0[YY])+sqr(x0[ZZ]))*massa[nn];
				Inertia[YY][YY]+=(sqr(x0[ZZ])+sqr(x0[XX]))*massa[nn];
				Inertia[ZZ][ZZ]+=(sqr(x0[XX])+sqr(x0[YY]))*massa[nn];
				Inertia[XX][YY]-=(x0[XX]*x0[YY])*massa[nn];
				Inertia[XX][ZZ]-=(x0[XX]*x0[ZZ])*massa[nn];
				Inertia[YY][ZZ]-=(x0[YY]*x0[ZZ])*massa[nn];
			}
		}
		Inertia[YY][XX]=Inertia[XX][YY];
		Inertia[ZZ][XX]=Inertia[XX][ZZ];
		Inertia[ZZ][YY]=Inertia[YY][ZZ];
		Jacob()(Inertia,Im,ei2,nrot);
		Im*=unit_nmm2;
		double fact=5.0/tmass/2.0;
		axis[XX]=fact*(Im[YY]+Im[ZZ]-Im[XX]);
		axis[YY]=fact*(Im[XX]+Im[ZZ]-Im[YY]);
		axis[ZZ]=fact*(Im[YY]+Im[XX]-Im[ZZ]);
		MyRg=(axis[XX]+axis[YY]+axis[ZZ])/5.0;
		Rg[o](MyRg,Im,Gm,axis);
	}

}
MAtoms::RelativeCoord::RelativeCoord(MAtoms *y,vectorvi & ind, vectorvi & mAtoms, string WhichAtoms){
	RelativeCoord(y,ind,mAtoms);
	if(WhichAtoms != "C") return;

}

MAtoms::RelativeCoord::RelativeCoord(MAtoms *y,vectorvi & ind, vectorvi & mAtoms){
	vector<Dvect> R_CM;
    vectorvd R_CM_m;
	Matrix CO(y->Mt.getCO());
	vector<double> & mass=y->mass;
	rvec * xa=y->xa;
	for(auto & it: ind){
	    Dvect cm=0.0;
	    double tmass=0.0;
	    vector<Dvect> cm_m;
		for(auto & p: it){
			Dvect cml_m=0.0;
			auto Mmass=0.0;
			for(auto q=0;q < mAtoms[p].size();q++){
				auto i=mAtoms[p][q];
	    		for(int o1=0;o1<DIM;o1++) {
	    			cm[o1]+=xa[i][o1]*mass[i];
	    			cml_m[o1]+=xa[i][o1]*mass[i];
	    		}
	    		tmass+=mass[i];
	    		Mmass+=mass[i];
			}
			cml_m/=Mmass;
			cm_m.push_back(cml_m);
		}
	    cm/=tmass;
	    R_CM.push_back(cm);
	    R_CM_m.push_back(cm_m);
	}
	for(auto i=0;i<R_CM.size();i++){
		Dvect cm=R_CM[i];
		for(auto j=0;j<R_CM_m[i].size();j++){
			R_CM_m[i][j]=R_CM_m[i][j]-R_CM[i];
		}
	}
	for(auto & it: R_CM)
		it=CO*it;
	for(auto & it: R_CM_m)
		for(auto & ip: it)
			ip=CO*ip;

	CM_c=R_CM;
	CM_m=R_CM_m;
}
void copy3x3AL(Matrix & x, alglib::real_2d_array & y){
	for(auto o=0;o<DIM;o++)
		for(auto p=0;p<DIM;p++)
			y[o][p]=x[o][p];
}
void copy3x3AL(alglib::real_2d_array & x, Matrix & y){
	for(auto o=0;o<DIM;o++)
		for(auto p=0;p<DIM;p++)
			y[o][p]=x[o][p];
}
Quaternions::Quaternion QfromMatrix(Matrix Q){
	const int x=0,y=1,z=2,p=3,dim=4;
	double t=Q[x][x]+Q[y][y]+Q[z][z];
	double r,s,ww,xx,yy,zz;
	if(t >= 0.0){
		r=sqrt(1+t);
		s=0.5/r;
		ww=0.5*r;
		xx=(Q[z][y]-Q[y][z])*s;
		yy=(Q[x][z]-Q[z][x])*s;
		zz=(Q[y][x]-Q[x][y])*s;
	} else {
		r=sqrt(fabs(1.0+Q[x][x]-Q[y][y]-Q[z][z]));
		s=0.5/r;
		ww=(Q[z][y]-Q[y][z])*s;
		xx=0.5*r;
		yy=(Q[x][y]+Q[y][x])*s;
		zz=(Q[z][x]+Q[x][z])*s;
	}
	Quaternions::Quaternion q={ww,xx,yy,zz};
	return q;

}
void MAtoms::CompRot(){

	CenterMass & R_cm=*R_cmx;
	vector<vector<int> > & mCluster=Perco->getCluster();
	vector<vector<int> > & mAtoms=Perco->getAtoms();
	alglib::real_2d_array A,U,VT,Out;
	alglib::real_1d_array W;
	A.setlength(3,3);
	Out.setlength(3,3);
	W.setlength(3);
	U.setlength(3,3);
	VT.setlength(3,3);
	CenterMass::setTime(time_c);
	if(!iCountRot) {
		vector<Quaternion> Q_t;

		auto myF=[&](int g){return Mdetg.find(atres[mAtoms[g][0]]) != string::npos;}; // find if detergent is present
		for(size_t o=0;o<mCluster.size();o++){
			int m=std::count_if(mCluster[o].begin(),mCluster[o].end(),myF);
			if(m){
				Ovec.push_back(vector<int>(m));
				Q_t.push_back(Quaternion{1.0,0.0,0.0,0.0});
				auto & Ot=Ovec.back();
				std::copy_if(mCluster[o].begin(),mCluster[o].end(),Ot.begin(),myF);
			}
		}
		RelativeCoord X(this,Ovec,mAtoms);
		CM_c0=X.getCM_c();
		CM_m0=X.getCM_m();
		CM_c1=CM_c0;
		CM_m1=CM_m0;
		try{ if(Q_t.size() != CM_c0.size()) throw " Micelles exist with no detergent. Cannot run ";}
		catch(const char * s){cout << s <<endl; exit(1);}
		R_cm.setup(Q_t.size()); R_cm(CM_c0,CM_m0,Q_t);

		Matrix CO(Mt.getCO());
		Matrix OC(Mt.getOC());
		R_cm.setCOs(CO,OC);
		iCountRot=1;

		R_cmx->RefCoord(x,atname,mCluster,mAtoms);
		return;
	}
	vector<Quaternion> Q_t(CM_c1.size());
	RelativeCoord X(this,Ovec,mAtoms);
	CM_c1=X.getCM_c();
	CM_m1=X.getCM_m();

	vector<Matrix> R(CM_m0.size());
	for(auto p=0;p<CM_m0.size();p++){
		R[p]=0.0;
		for(auto q=0;q<CM_m0[p].size();q++){
			R[p]+=CM_m1[p][q]%CM_m0[p][q];
		}
		R[p]/=static_cast<double>(CM_m0[p].size());
		copy3x3AL(R[p],A);
		alglib::rmatrixsvd(A,3,3,2,2,2,W,U,VT);
		alglib::rmatrixgemm(3, 3, 3, 1, U, 0, 0, 0, VT, 0, 0, 0, 0.0, Out, 0, 0);
		copy3x3AL(Out,R[p]);
		Quaternion q=QfromMatrix(R[p]);
		Q_t[p]=q;
	}

	R_cm(CM_c1,CM_m1,Q_t);
	Matrix CO(Mt.getCO());
	Matrix OC(Mt.getOC());
	R_cm.setCOs(CO,OC);
	CenterMass::setTime(time_c);
	CM_c0=CM_c1;
	CM_m0=CM_m1;
	R_cmx->RefCoord(x,atname,mCluster,mAtoms);

	/*
		WignerDMatrix Wigner(Q_t[p]);
		vector<complex<double> > Wig;

		for(auto i=-1;i<=1;i++)
			for(auto j=-1;j<=1;j++)
				Wig.push_back(Wigner(1,i,j));
		for(auto i=-2;i<=2;i++)
			Wig.push_back(Wigner(2,i,0));

	 */
}
void MAtoms::CompCM(){
	CenterMass & R_cm=*R_cmx;
	vector<vector<int> > & mCluster=Perco->getCluster();
	vector<vector<int> > & mAtoms=Perco->getAtoms();
	vector<Dvect> R_CM=vector<Dvect>(mCluster.size());
	CenterMass::setTime(time_c);
	for(size_t o=0;o<mCluster.size();o++){
		Dvect cm=0.0;
		double tmass=0.0;
		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];
			for(size_t q=0;q<mAtoms[n].size();q++){
				int i=mAtoms[n][q];
				for(int o1=0;o1<DIM;o1++) cm[o1]+=xa[i][o1]*mass[i];
				tmass+=mass[i];
			}
		}
		cm/=tmass;
		R_CM[o]=cm;
	}
	R_cm(R_CM);
	double a=Mt.getCO()[XX][XX];
	double b=Mt.getCO()[YY][YY];
	double c=Mt.getCO()[ZZ][ZZ];
	R_cm.setAxis(a,b,c);
  //  copy(R_CM.begin(),R_CM.end(),back_inserter(R_cm));
}

void MAtoms::Gyro(){
	vector<vector<int> > mCluster=Perco->getCluster();
	vector<vector<int> > mAtoms=Perco->getAtoms();
	vector<Gyration> Rg=vector<Gyration>(mCluster.size());
	Rg_i.clear();
	Gyration::setTime(time_c);
	CalcGyro(mass,Rg);
	copy(Rg.begin(),Rg.end(),back_inserter(Rg_i));
	CalcGyro(massNCH,Rg);
	copy(Rg.begin(),Rg.end(),back_inserter(Rg_i));

	Rg_count++;
}
Dvect MAtoms::CM(){
	Dvect cm=0.0;
	double tmass=0.0;
	for(unsigned int n=0;n<ResIndx0->size();n++){
		int nn=(*ResIndx0)[n];
		for(int o=0;o<DIM;o++) cm[o]+=x[nn][o]*mass[nn];
		tmass+=mass[nn];
	}
	cm/=tmass;
	return cm;
}
Dvect MAtoms::CMCenter(){
	Dvect cm=0.0;
	double tmass=0.0;
	for(unsigned int n=0;n < MyRes.size();n++){
		for(unsigned int m=0;m < MyRes[n].size();m++){
			int nn=MyRes[n][m];
			for(int o=0;o<DIM;o++) cm[o]+=x[nn][o]*mass[nn];
			tmass+=mass[nn];
		}
	}
	cm/=tmass;
	return cm;
}
Matrix MAtoms::InertiaRotCenter(){
	Dvect x0,cm;
	cm=CMCenter();

	double unit_nmm2=1.0/(unit_nm*unit_nm);
	Matrix Inertia=0.0,ei,eiT;
	Dvect Im;
	int nrot;

	for(unsigned int n=0;n < MyRes.size();n++){
		for(unsigned int m=0;m < MyRes[n].size();m++){
			int nn=MyRes[n][m];
			for(int o=0;o<DIM;o++) x0[o]=x[nn][o]-cm[o];
			Inertia[XX][XX]+=(sqr(x0[YY])+sqr(x0[ZZ]))*mass[nn];
			Inertia[YY][YY]+=(sqr(x0[ZZ])+sqr(x0[XX]))*mass[nn];
			Inertia[ZZ][ZZ]+=(sqr(x0[XX])+sqr(x0[YY]))*mass[nn];
			Inertia[XX][YY]-=(x0[XX]*x0[YY])*mass[nn];
			Inertia[XX][ZZ]-=(x0[XX]*x0[ZZ])*mass[nn];
			Inertia[YY][ZZ]-=(x0[YY]*x0[ZZ])*mass[nn];
		}
	}
	Inertia[YY][XX]=Inertia[XX][YY];
	Inertia[ZZ][XX]=Inertia[XX][ZZ];
	Inertia[ZZ][YY]=Inertia[YY][ZZ];
	Jacob()(Inertia,Im,ei,nrot);
	eiT=ei.Transpose();
	return eiT;
}

Matrix MAtoms::InertiaRot(){
	Dvect cm=0.0,x0;
	double tmass=0.0;
	for(unsigned int n=0;n<ResIndx0->size();n++){
		int nn=(*ResIndx0)[n];
		for(int o=0;o<DIM;o++) cm[o]+=x[nn][o]*mass[nn];
		tmass+=mass[nn];
	}
	cm/=tmass;

	double unit_nmm2=1.0/(unit_nm*unit_nm);
	Matrix Inertia=0.0,ei,eiT;
	Dvect Im;
	int nrot;
	for(unsigned int n=0;n<ResIndx0->size();n++){
		int nn=(*ResIndx0)[n];
		for(int o=0;o<DIM;o++) x0[o]=x[nn][o]-cm[o];

		Inertia[XX][XX]+=(sqr(x0[YY])+sqr(x0[ZZ]))*mass[nn];
		Inertia[YY][YY]+=(sqr(x0[ZZ])+sqr(x0[XX]))*mass[nn];
		Inertia[ZZ][ZZ]+=(sqr(x0[XX])+sqr(x0[YY]))*mass[nn];
		Inertia[XX][YY]-=(x0[XX]*x0[YY])*mass[nn];
		Inertia[XX][ZZ]-=(x0[XX]*x0[ZZ])*mass[nn];
		Inertia[YY][ZZ]-=(x0[YY]*x0[ZZ])*mass[nn];
	}
	Inertia[YY][XX]=Inertia[XX][YY];
	Inertia[ZZ][XX]=Inertia[XX][ZZ];
	Inertia[ZZ][YY]=Inertia[YY][ZZ];
	Jacob()(Inertia,Im,ei,nrot);
	eiT=ei.Transpose();
	return eiT;
}
string MAtoms::gTotResList(){
	stringstream ss;
	for(vector<string>::iterator it=ResList0->begin();it != ResList0->end();it++)
		ss<< *it + " " ;
	return ss.str();
}
void MAtoms::InitSelection(vector<string> & y, TopolMic & MyTop){
	try{
		for(size_t o=0;o<y.size();o++){
			if(!MyTop.CheckResidue(y[o])) throw string(" Reference residue " + y[o] + " does not exist in the system. Abort");
		}
	} catch(const string & s){
		cout << s << endl;
		exit(1);
	}
	typedef map<string,vector<vector<int> > >  ResMap;
	typedef ResMap::iterator ResIter;

	ResMap & ResDef=MyTop.getDef();
	for(size_t o=0;o<y.size();o++){
		vector<vector<int> > Res=ResDef[y[o]];
		for(size_t p=0;p<Res.size();p++)
			SelRes.push_back(Res[p]);
	}
}
void MAtoms::Reconstruct(const string & y, TopolMic & MyTop){
	try{
		if(!MyTop.CheckResidue(y)) throw string(" Reference residue " + y + " does not exist in the system. Abort");
	} catch(const string & s){
		cout << s << endl;
		exit(1);
	}
	if(!calls++) {PBCvect dummy;}
	vector<Dvect> nxyz=PBCvect::getVec();
	typedef map<string,vector<vector<int> > >  ResMap;
	typedef ResMap::iterator ResIter;

	ResMap & ResDef=MyTop.getDef();
	vector<vector<int> > & Res=ResDef[y];
	MyRes=ResDef[y];
	Matrix co=Mt.getCO();
	Matrix oc=Mt.getOC();
	double v[DIM];
	v[XX]=co[XX][XX];v[YY]=co[YY][YY];v[ZZ]=co[ZZ][ZZ];
	double MinCO=*min_element(v,v+3)*0.5;  // Half the smallest axis is the cutoff distance for atom shifts
	Dvect Xref;

	// Apply boundary condition on every molecule as the default boundaries n .xtc files are atomic

	for(ResIter it=ResDef.begin(); it != ResDef.end();it++){
		vector<vector<int> > & Res0=it->second;
		for(unsigned int o=0;o<Res0.size();o++){
			Xref=xa[Res0[o][0]];
			for(unsigned int p=1;p<Res0[o].size();p++) {
				Dvect x1=Xref-xa[Res0[o][p]];
				Dvect xc1=co*x1;
				if(xc1.Norm() > MinCO) {
					Dvect tmp=*min_element(nxyz.begin(),nxyz.end(),CellComp(x1,co));
					for(int q=0; q < DIM;q++) xa[Res0[o][p]][q]=xa[Res0[o][p]][q]-tmp[q];
				}
			}
		}
	}
	// > Take the reference type molecules and eliminate boundary conditions --> begin

	// >>  Compute center of mass for each molecule


	vector<Dvect> xaa;
	for(unsigned int o=0; o< Res.size();o++){
		Dvect cm=0.0;
		for(unsigned int p=0; p< Res[o].size();p++){
			int pp=Res[o][p];
			cm+=xa[pp];
		}
		cm/=static_cast<double>(Res[o].size());
		xaa.push_back(cm);
	}


	// >> Obtain the center of mass for the cluster of the reference molecules

	Xref=xaa[0];
	for(unsigned n=1;n < xaa.size();n++){
		Dvect xb=xaa[n];
		Dvect tmp=*min_element(nxyz.begin(),nxyz.end(),CellComp(Xref,xb,co));
		xaa[n]=xaa[n]-tmp;
	}

	Dvect cma=0.0;
	for(unsigned n=0;n < xaa.size();n++){
		cma+=xaa[n];
	}
	cma/=static_cast<double> (xaa.size());
	// >>Shift all system by cma

	for(int n=0;n<nr;n++)
		for(int m=0;m<DIM;m++)
			xa[n][m]=xa[n][m]-cma[m];


	//> Apply molecular boundary conditions on the entire system
	for(ResIter it=ResDef.begin(); it != ResDef.end();it++){
		vector<vector<int> > & Res0=it->second;
		for(unsigned int o=0;o<Res0.size();o++){
			Dvect xcm=0.0;
			for(unsigned int p=0;p< Res0[o].size();p++) {
				xcm+=xa[Res0[o][p]];
			}

			xcm/=static_cast<double>(Res0[o].size());

			for(int m=0;m<DIM;m++) xcm[m]=rint(xcm[m]);
			for(unsigned int p=0;p<Res0[o].size();p++)
				for(int m=0;m<DIM;m++)
					xa[Res0[o][p]][m]=xa[Res0[o][p]][m]-xcm[m];

		}
	}

		//> Obtain new Cartesian coordinates from reduced xa's
	for(int i=0;i<nr;i++){
		for(int o=0;o<DIM;o++){
			xa[i][o]=xa[i][o]+0.5;
			x[i][o]=Mt.getCO()[o][XX]*xa[i][XX]+Mt.getCO()[o][YY]*xa[i][YY]+Mt.getCO()[o][ZZ]*xa[i][ZZ];
		}
	}
}
void MAtoms::Reconstruct(){
	if(!calls++) {PBCvect dummy;}
	vector<Dvect> nxyz=PBCvect::getVec();

	vector<vector<int> > mCluster=Perco->getCluster();
	vector<vector<int> > mAtoms=Perco->getAtoms();
	Matrix co=Mt.getCO();
	Matrix oc=Mt.getOC();
	double v[DIM];
	v[XX]=co[XX][XX];v[YY]=co[YY][YY];v[ZZ]=co[ZZ][ZZ];
	double MinCO=*min_element(v,v+3)*0.5;  // Half the smallest axis is the cutoff distance for atom shifts
	Dvect Xref;

	// Loop on the Clusters


	for(size_t o=0;o<mAtoms.size();o++){
		Xref=xa[mAtoms[o][0]];
		for(size_t p=0;p<mAtoms[o].size();p++){
			int n=mAtoms[o][p];
			Dvect x1=Xref-xa[n];
			Dvect xp=xa[n];
			Dvect xc1=co*x1;
			if(xc1.Norm() > MinCO) {
				Dvect tmp=*min_element(nxyz.begin(),nxyz.end(),CellComp(x1,co));
				for(int q=0; q < DIM;q++) xa[n][q]=xa[n][q]-tmp[q];
			}

		}
	}

	for(size_t o=0;o<mCluster.size();o++){
		Xref=0.0;

		vector<Dvect> xcm=vector<Dvect>(mCluster[o].size(),0.0);
		Dvect xcmC=0.0;

		// Compute molecule center of mass xcm

		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];

			// Apply boundary condition on every molecule as the default boundaries n .xtc files are atomic
			for(size_t i=0;i<mAtoms[n].size();i++){
				Dvect xx=xa[mAtoms[n][i]];
				xcm[p]+=xx;
			}
			xcm[p]/=static_cast<double>(mAtoms[n].size());
		}

		Xref=xcm[0];  // Reference to the first molecule of the cluster



		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];
			Dvect x1=Xref-xcm[p];
			Dvect xcmp=xcm[p];
			Dvect xc1=co*x1;

			if(xc1.Norm() > MinCO) {
				Dvect tmp=*min_element(nxyz.begin(),nxyz.end(),CellComp(x1,co));
				for(int q=0; q < DIM;q++) xcm[p][q]=xcm[p][q]-tmp[q];

				for(size_t i=0;i<mAtoms[n].size();i++){
					int ia=mAtoms[n][i];
					for(int q=0; q < DIM;q++) xa[ia][q]=xa[ia][q]-tmp[q];
				}
			}

		}
		// Compute the center of mass of the cluster

		for(size_t p=0;p<mCluster[o].size();p++){
			xcmC+=xcm[p];
		}
		xcmC/=static_cast<double>(mCluster[o].size());

		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];
			for(size_t i=0;i<mAtoms[n].size();i++){
				int ia=mAtoms[n][i];
				xa[ia][XX]=xa[ia][XX]-rint(xcmC[XX]);
				xa[ia][YY]=xa[ia][YY]-rint(xcmC[YY]);
				xa[ia][ZZ]=xa[ia][ZZ]-rint(xcmC[ZZ]);
			}

		}

	}

	//> Obtain new Cartesian coordinates from reduced xa's

	for(size_t i1=0;i1<mCluster.size();i1++){
		for(size_t i2=0;i2<mCluster[i1].size();i2++){
			int n=mCluster[i1][i2];
			for(size_t i=0;i<mAtoms[n].size();i++){
				int ia=mAtoms[n][i];
				for(int o=0;o<DIM;o++){
					xa[ia][o]=xa[ia][o]+0.5;
					x[ia][o]=Mt.getCO()[o][XX]*xa[ia][XX]+Mt.getCO()[o][YY]*xa[ia][YY]+Mt.getCO()[o][ZZ]*xa[ia][ZZ];
				}
			}
		}
	}
}

MAtoms::~MAtoms() {
	delete Perco;
	delete R_cmx;
	// TODO Auto-generated destructor stub
}

ostream & operator<<(ostream & fout, MAtoms & y){
	y.cPrint(fout);
	return fout;
}

void MAtoms::cPrint(ostream & fout){
	vector<Dvect> xc=vector<Dvect>(this->nr);

	Matrix Rot0=InertiaRotCenter();
	Dvect cm=CMCenter();

	for(int n=0;n<nr;n++){
//		Dvect xb=x[n];
		xc[n]=x[n];
//		xb=xb-cm;
//		xc[n]=Rot0*xb;
	}


	stringstream ss;
	ss<<time_c;
	fout << string("REMARK    Generated by trjMicelles ") <<endl;
	fout << string("REMARK    SIMULATION TIME = "+ ss.str()) <<endl;

	vector<double> Par=Mt.getParas();
	ss.str(string());
	ss<< setw(9) << setprecision(3) << fixed << Par[0]*10.0;
	ss<< setw(9) << setprecision(3) << fixed << Par[1]*10.0 ;
	ss<< setw(9) << setprecision(3) << fixed << Par[2]*10.0 ;
	ss<< setw(7) << setprecision(2) << fixed << Par[3] ;
	ss<< setw(7) << setprecision(2) << fixed << Par[4] ;
	ss<< setw(7) << setprecision(2) << fixed << Par[5] ;
	fout<< string("CRYST1"+ ss.str() + " P 1           1") <<endl;

	map<int,vector<string> >MapRes;
	for(unsigned int n=0;n<ResIndx0->size();n++){
		string tmp;
		int o=(*ResIndx0)[n];
		ss.str(string());
		ss <<std::setw(5) << std::left << PDB[o].resn;
		string tmp2=PDB[o].first;
		tmp2.erase(26,1);
		tmp=tmp2.replace(21,5,ss.str());
		ss.str(string());
		ss << std::setw(5) << std::right << PDB[o].res+1;
		tmp.replace(21,5,ss.str());
		tmp.append(" 1.00");
		tmp.append(PDB[o].last);
		ss.str(string());
		ss << fixed << setw(8) << setprecision(3) << right << xc[o][XX]/unit_nm;;
		ss << fixed << setw(8) << setprecision(3) << right << xc[o][YY]/unit_nm;;
		ss << fixed << setw(8) << setprecision(3) << right << xc[o][ZZ]/unit_nm;;
		tmp.replace(30,24,ss.str());
		MapRes[PDB[o].res+1].push_back(tmp);
	}

	map<int,vector<string> >::iterator itt=MapRes.begin();
	for(;itt != MapRes.end(); ++itt){
		vector<string> & tmp=itt->second;
		std::sort(tmp.begin(),tmp.end(),op_sortAtnoPDB());
		for(size_t o=0; o < tmp.size(); o++){
			fout <<tmp[o] <<endl;
		}
	}
	fout << string("TER") <<endl;
	fout << string("ENDMFL") <<endl;
}
void MAtoms::fPrint(Fstream & fout){
	FstreamC * foutC=dynamic_cast<FstreamC *> (&fout);
	t_fileio *xd = foutC->getfin();
	gmx_bool bOK=bOK_c;
	int step=step_c;
	float prec=prec_c;
	float time=time_c;
	matrix box;

	vector<Dvect> x0=vector<Dvect>(ResIndx0->size());
	for(size_t i=0;i < x0.size();i++){
		int ia=ResIndx0->at(i);
		x0.push_back(x[ia]);
	}
	int    natos=static_cast<int> (x0.size());
	rvec * x=new rvec [natos];
	for(unsigned int m=0;m<x0.size();m++)
		for(int o=0;o<DIM;o++) x[m][o]=x0[m][o];

	for(int i=0;i<DIM;i++)
			for(int j=0;j<DIM;j++)
				box[i][j]=Mt.getCO()[i][j];
	try{
		if(!write_xtc(xd,natos,step,time,box,x,prec)) throw string("Cannot write next frame. Abort");
	}catch(const string & s){
		cout << s << endl;
		exit(1);
	}

}
Fstream & operator<<(Fstream & fout, MAtoms & y){
	y.fPrint(fout);
	return fout;
}


